%load policy functions from risk-free debt model
cd ..
cd ..
cd ..
cd 'riskfree_debt/defaus/output'

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');

nb2=3501;

b   =load('bgrid2.txt');
b   = b/.204;

for i=1:na;
    for j=1:nb2;
        vconrf(i,j)  =vcon0((i-1)*nb2 + j);
        
        prf(i,j)  =p0((i-1)*nb2 + j);
        hrf(i,j)  =h0((i-1)*nb2 + j);
        ibprf(i,j)=ibp0((i-1)*nb2 + j);
        ctrf(i,j) =ct0((i-1)*nb2 + j);
        carf(i,j) =ca0((i-1)*nb2 + j);
        gnrf(i,j) =gn0((i-1)*nb2 + j);       
        taurf(i,j)=tau0((i-1)*nb2 + j);       
        wrf(i,j)  =w0((i-1)*nb2 + j);   
        defrf(i,j)=def0((i-1)*nb2 + j);   
        qrf(i,j)  =q0((i-1)*nb2 + j);   
        recrf(i,j)=wrf(i,j)*hrf(i,j)*taurf(i,j);
    end;
end;

bprf = b(ibprf);

for i=1:na;
    vautrf(i)  =vaut0(i);
    
    pautrf(i)  =paut0(i);
    hautrf(i)  =haut0(i);
    ctautrf(i) =ctaut0(i);
    gnautrf(i) =gnaut0(i);       
    tauautrf(i)=tauaut0(i);       
    wautrf(i)  =waut0(i);   
    recautrf(i)=waut(i)*haut(i)*tauaut(i);
	wwautrf(i) = alpha*paut(i).*haut(i).^(alpha-1); 
    
end;

cnrf       = hrf.^alpha-gnrf;
cnautrf    = hautrf.^alpha-gnautrf;
constrf    = (omegac.*ctrf.^(-mu)+(1-omegac).*cnrf.^(-mu)).^(-1/mu);
constautrf = (omegac.*ctautrf.^(-mu)+(1-omegac).*cnautrf.^(-mu)).^(-1/mu);

for i=1:na;
    temp1 = find(ctrf(i,:) > 0);
    temp2 = find(cnrf(i,:) > 0);
    temp3 = find(prf(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    vrf(i,1:imin(i)-1)  = NaN; 
    prf(i,1:imin(i)-1)  = NaN; 
    bprf(i,1:imin(i)-1) = NaN; 
    ibprf(i,1:imin(i)-1)= NaN; 
    ctErf(i,1:imin(i)-1) = NaN; 
    ctUrf(i,1:imin(i)-1) = NaN; 
    cnErf(i,1:imin(i)-1) = NaN; 
    cnUrf(i,1:imin(i)-1) = NaN; 
    constErf(i,1:imin(i)-1) = NaN;     
    constUrf(i,1:imin(i)-1) = NaN;     
    gnrf(i,1:imin(i)-1) = NaN; 
    hrf(i,1:imin(i)-1)  = NaN; 
    carf(i,1:imin(i)-1)  = NaN; 
    taurf(i,1:imin(i)-1)= NaN; 
	recrf(i,1:imin(i)-1)= NaN;
	wwrf(i,1:imin(i)-1)= NaN;
	qrf(i,1:imin(i)-1)  = NaN;
  
end;

